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EXECUTIVE  SUMMARY 


Introduction 

Detailed,  quantitative  interpretations  of  images  require  data  about  the  angular  distribution  of  light 
reflected  from  a  shadowed  sea.  This  report  demonstrates  the  solution  of  the  two-dimensional 
radiative  transfer  problem  for  the  sea  shadowed  by  an  opaque  body  such  as  a  pier  (See  Figure  1). 
The  solution  to  the  problem  requires  two  different  theoretical  approaches.  The  solution  for  the 
depth  dependence  of  light  radiance  is  obtained  using  the  theory  of  Green’s  functions  (Vladimirov, 
1971).  The  solution  to  the  problem  in  horizontal  direction  is  obtained  with  the  one-dimensional 
Fourier  transform  method  (Morse  and  Feshbach,  1953). 

The  algorithm  used  in  this  report  is  based  on  the  analytical  solution  to  the  spatially  two-dimensional 
(Cartesian  coordinates  x  and  z)  scalar  radiative  transfer  equation  in  seawater.  The  Ox  axis  starts  on 
toe  sea  surface  in  toe  center  of  toe  opaque  shadow  (pier)  and  increases  towards  toe  sun  in  toe 
direction  orthogonal  to  toe  axis  of  the  pier.  The  Oz  axis  increases  down  to  toe  sea  bottom. 

Because  toe  input  parameters  did  not  include  toe  scattering  phase  function,  toe  algorithm  uses 
experimental  results  published  in  literature.  The  scattering  phase  function  is  estimated  through  toe 
regression  relationships  which  are  derived  by  processing  aU  fifteen  experimental  marine  Petzold 
phase  functions  (Haltrin,  1997).  This  regression  relationship  couples  toe  scattering  phase  function 
with  toe  scattering  coefficient  and  toe  single-scattering  albedo. 

The  future  possible  development  of  this  problem  may  consists  of  toe  following  enhancements  to 
toe  theory:  (a)  consideration  of  toe  reflection  by  toe  sea  bottom;  (b)  consideration  of  toe 
atmospheric  scattering  effects  (Haltrin,  1998b,  1998c)  and  toe  elevation  of  toe  pier;  (c) 
consideration  of  toe  sea  surface  roughness;  (d)  generalization  to  toe  non-infinite  pier  (a  3D-problem 
in  space);  (e)  consideration  of  toe  vertical  optical  stratification  of  seawater;  and  (f)  enhancement  of 
toe  part  of  the  theory  related  to  toe  higher-order  scattering  effects. 


Findings  and  Recommendations 

The  final  equation  obtained  in  this  report  expresses  toe  upwelling  sea  radiance  distribution  as  a 
function  of  toe  sun  zenith  angle,  azimuth  and  zenith  angles  of  viewing,  toe  width  of  toe  pier  and 
toe  inherent  optical  properties  of  toe  water  (absorption  and  attenuation  coefficients  and  toe 
scattering  phase  function).  Results  for  toe  actual  radiance  distribution  of  light  are  obtained 
numerically  with  toe  program  RADPIER  (see  Appendix)  using  as  input  toe  in  situ  measured  optical 
properties  of  seawater.  The  scattering  phase  function  is  estimated  through  toe  regression  equations 
that  are  based  on  aU  fifteen  Petzold  scattering  phase  functions.  The  empirically  derived  scattering 
phase  function  depends  on  toe  scattering  coefficient,  toe  single-scattering  albedo,  and  toe  scattering 
angle.  The  resulting  radiance  distributions  are  discussed  and  illustrated. 
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INTRODUCTION 

Our  goal  is  formulated  as  follows:  Given  a  downwelhng  radiance  distribution  from  the  sun 
consistent  with  a  long  rectangular  shadow,  to  model  the  water-leading  radiance  as  a  function  of 
azimuth  and  zenith  angles  of  viewing  and  distance  in  the  direction  normal  to  the  long  dimension  of 
the  rectangle.  The  derived  approach  should  be  vahd  for  an  arbitrary  solar  elevation. 

Several  serious  attempts  have  been  made  to  solve  two-  and  three-dimensional  radiative  transfer 
problem  (Vladimirov,  1961;  Williams,  1967;  Kochetkov,  1968;  Kaper,  1969;  Case  and  Hazeltine, 
1970,  Plateau  and  Stephens,  1988).  Unfortunately  none  of  these  approaches  are  simple  emd 
convenient  enough  to  produce  practically  usable  equations  applicable  to  solve  our  problem. 


FORMAL  METHODS 

We  start  from  the  scalar  radiative  transfer  theory  (Chandrasekhar,  1960).  Our  problem  is  a  two 
dimensional  in  space  radiative  transfer  problem.  TTie  origin  of  the  Cartesian  coordinates  is  placed 
on  the  sea  surface  just  below  the  center  of  an  infinite  pier.  For  simplicity  we  place  the  pier 
horizontally  along  the  Oy  axis,  and  the  solar  plane  is  perpendicular  to  the  Oy  axis.  The  Ox  axis  is 
directed  horizontally  to  Ae  sun  and  ties  in  the  solar  plane,  and  Oz  axis  is  directed  vertically  to  the 
bottom. 

We  use  the  Green’s  function  method  (Vladimirov,  1971)  to  solve  the  problem  in  a  z-direction,  and 
the  Fourier  transform  method  to  solve  an  x-dependent  part  (Morse  and  Feshbach,  1953). 

The  angular  distribution  of  radiance  is  split  into  three  components:  unscattered,  single-scattered  and 
multiple  scattered  components.  Exact  solutions  for  the  first  two  components  are  obtained. 


HOMOGENEOUS  SHALLOW  SEA  ILLUMINATED  BY  THE  SUNLIGHT 

We  assume  that  the  reader  is  familiar  with  the  standard  terminology  of  hydrologic  optics  (see,  for 
example,  Jerlov,  1986;  Shifrin  1988,  Mobley,  1994). 

Let  us  introduce  the  following  notations:  is  irradiance  by  the  sun  on  the  sea  surface, 

is  irradiance  by  the  sun  on  the  surface  normal  to  the  solar  rays,  fj.^  =  cosZ^  =  sinh^  is 
the  cosine  of  the  solar  zenith  angle  Z^,  hg  =  90°  -  Z^  is  a  solar  elevation  angle. 

The  solar  radiance  above  the  sea  surface  is  expressed  by: 

Ll  =  E°6in-nl)5((p),  (1) 

where  5(x)  is  the  Dirac’s  delta-function  (Morse  and  Feshbach,  1953) 
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is  the  cosine  of  the  light  penetration  angle,  is  the  seawater  refraction  coefficient. 

If  r/  =  i-/?i(z,)  is  a  transmission  of  light  by  sea  surface,  where  is  a  Fresnel  reflection 
coefficient  of  light  falling  from  above,  than  the  light  radiance  below  the  sea  surface  is: 

=  £!.  =  £^r/,  (2) 

here  is  the  solar  irradiance  on  the  surface  normal  to  the  rays  just  below  the  sea  surface. 

A  one-dimensional  scalar  radiative  transfer  equation  (Chandrasekhar,  1960)  for  the  total  light 
radiance  can  be  written  as: 


d 

u  —  +  c 
dz 


L(z,H,(p)  =  —  jd(p'jpicosy)^z,^l\<p')dn', 


(3) 


here  c  =  a  +  b  is  the  attenuation  coefficient,  a  is  the  absorption  coefficient  and  b  is  the  scattering 
coefficient  of  seawater,  fi  =  cos0,  the  direction  of  Ught  propagation  is  determined  by  the  zenith 
and  azimuth  angles  6  and  <p,  p(cosy)  is  the  scattering  phase  function,  /  is  the  scattering  angle 
determined  from  the  following  formula: 

cosY  =  ^l^l'  +  ^|l-p^^|\-  cos(<p  -  (p') .  (4) 

Let  us  represent  the  total  radiance  of  light  inside  the  sea  L  as  the  direct  light  radiance  1^,  and  the 
sum  of  radiances  that  represent  order  of  scattering: 


£  =  (5) 

n=0 

By  substituting  Eq.  (5)  into  Eq.  (3)  we  obtain  the  following  equations  for  radiance  components  of 
different  orders  of  scattering  denoted  by  a  subscript: 

(  d  \ 

p  —  -^c\L^(z,H,(p)  =  0,  (6) 

\  dz  ) 


DL„{z,p,<p)  =  bSL„_^iz,p,(p),  or  L„{z,H,q>)  =  bD^'S L„^i(z,H,(p)  =  bf  L„_,iz,p,<p).  (7) 
The  operators  D  and  5  are  defined  as  follows: 

DF(z,p,(p)  =  +  c^F(z,H,<p), 

.  2n  \ 

S  F(z,p,<p)  =  —  \d(p'\ p(cosY)F(z,n\(p')dp', 


-1 


(8) 
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and  the  following  relationships  are  vaUd: 

^  .  At  A.  AA 

Z)‘‘=G,  T  =  D-^S  =  GS. 


(8a) 


_  A 

The  operator  G  in  Eq.  (8a)  is  a  Green’s  function  of  Eq.  (3).  This  operator  is  defined  according  to 
the  following  formula  (Vladimirov,  1971): 


•foo 

GF{z)^jG{z-z)Fiz')dz', 


(9) 


Here  H{z)  is  the  Heavyside’s  or  step  function  {H{z)  =  1,  z  ^  0;  H{z)  =  0,  z  <  0)  .(Morse  and 
Feshbach,  1953).  The  total  scattering  operator  is: 


^  2n  \ 

fF{z,^L,(p)  =  —  \d(p'\d^'pi,cosy)\dz'G{z-z')F{z\F'^(p')  •  (10) 


0 

Now,  the  total  radiance  is  expressed  as  a  following  sum: 


L  =  ^b’‘L„,  L„=F'rL,,  «>1,  L(z,H,(p)  =  %b'‘rL,{z,H,(p). 


(11) 


B  =  0 


n=0 


Let  us  proceed  to  solve  these  equations. 


SOLUTIONS  FOR  THE  RADIANCE 

The  solution  for  direct  (unscattered)  radiance  is  straightforwardly  obtained  from  the  Eq.  (6) 
with  the  following  boundary  condition  on  the  sea  surface,  L^{0,^,(p)  =  : 


L^{z,li,(p)  =  El  -  fis)S{<p) . 

Equation  (7)  for  the  single-scattered  radiance  can  be  written  as: 


(12) 


^  'Ik  1  -roo 

L,(z,m,<P)  ^  —  j d(p'\ dll' picosy)  J dz'G{z - z')hiz',li\¥)  ■  (13) 


0  -1 


By  substituting  Eq.  (12)  into  Eq.  (13),  we  have  the  following  solution  for  the  single-scattered 
radiance  distribution: 


mz,H,(p)=  -  [\l/iiz^li)H(ii)  +  \i/2iz,li)Hi-ii)], 

Ancfi 


(14) 
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g-CZlfls  _g-czlli 

¥iiz,M)  = - J - 1 - ,  =  W\(z,H)l^g  =  cz,  (15) 


W2iz,H) 


_1_  ]_ 


(16) 


For  the  required  precision  about  15-20%  it  is  sufficient  (see  Shuleykin,  1968)  to  calculate  the  next 
term  in  expansion  given  by  Eq.  (11).  The  explicit  solution  for  the  next  term  can  be  obtained  from 
the  following  equation: 


*  In  \  -H» 

l^{z,pi,(p)  =  —  \d(p'^dii'p{cosY)^dzG{z-z)mz,n',(p') .  (17) 

^^0-1  -<» 

with  L,  given  by  Eqs.  (14)-(16)  and  G(z)  given  by  Eq.  (9).  The  right  boundary  conditions  has  no 
influence  on  upwelling  radiance  distribution,  so  we  do  not  discuss  them  here. 


SOLUTIONS  FOR  THE  UPWELLING  RADIANCE  NEAR  THE  SEA  SURFACE 


Taking  integrals  in  (17)  and  summing  all  required  terms  at  z  =  0  and  /t  <  0,  we  have  the  following 
solution  for  the  case  of  homogeneous  illumination: 


here 


Li0,^,(p)  = 


El^sb 


47ic(^,  +  \n\)l 


p(cosy,)  +  —  \i/^iH,<p) 


,  M<0, 


cosYs  =  HHs  +  cosip , 


y/p{H,(p)  =  l/t|  f  J X  (p,p\<p)dp\ 

Xp{p,n',(p)  =  —  Jp(cosy)p(cosy;)d<p', 

ZTT  Q 

cosy;  =  fi%  +  cosip' . 


(18) 

(19) 

(20) 

(21) 

(22) 
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HOMOGENEOUS  SEA  ILLUMINATED  BY  SUNLIGHT  WITH  SHADOW 


Let  us  formulate  a  two-dimensional  problem  that  takes  into  account  an  unhomogeneous  over 
horizontal  axis  Ox  illumination.  The  two-dimensional  radiative  transfer  equation  for  the  total 
angular  radiance  distribution  is: 


f  I 

dz 


d  r.  T  d 

jU— +  ^l-/x  COS(p—  +  C 


L{x,  z,l2,(p)  =  ^jd(p'j  jp(cos  y)L{x,z,n',(p')diii',  (23) 


271  1 


0  -1 


Let  us  define  a  one-dimensional  Fourier  transform  by  the  following  equations: 

-foo  +00 

=  nx)  =  ^lf,e‘‘dk. 


(24) 


After  taking  a  Fourier  transform  Eq.  (23)  transfers  to  the  one-dimensional  equation  for  the  Fourier 
amplitude  L^(z,lii,(p): 


here 


2n  1 


(U —  +  C  J  4  (z,  /t,  <p)  =  —  J  J  p(cos  7)  4  (z,  n',  <p')  d/i'. 


47t 


-1 


C  =  c(l  +  ikT),  T  =  ^jl-^  cos  (pjc. 


(25) 

(26) 


Equation  (25)  may  be  formally  obtained  from  Eq.(3)  by  replacing  the  real  extinction  coefficient  c 
by  the  complex  value  c  given  by  Eq.  (26).  It  means  that  the  solution  of  Eq.  (24)  can  be  obtained 
from  the  solution  (18)  by  replacing  ^  x-dependent  values  by  their  Fourier  transforms  and  aU 
instances  of  the  extinction  coefficient  c  by  the  c.  Taking  this  into  account,  we  have  the  following 
solution  to  Eq.  (25): 


4(0,/i,<p)  = 


^wk  f^S  ^ 
47t  c{iii,  +  \n\) 


p(cosys)  +  —Wpi^i,(P) 

2c 


(27) 


Now  we  only  need  to  calculate  a  Fourier  transform  of  the  x-dependent  radicuice  distribution 


+0O 


(28) 


In  our  case  the  solar  illumination  incorporates  a  shadowing  at  -w  <  x  <  w ,  where  w  is  a  half 
width  of  a  shadow.  The  angular-space  distribution  of  the  hght  just  under  the  sea  surface  is: 

Llip,(p)  =  E^(x)d(^-Hs)Si(p)>  (29) 

E,,(x)  =  El[l-  H(w  -  |a:|)]  =  {l  -  0.5[sign(H'  -  x)  -h  sign(w  -t-  x)]} , 

here  sign(jc)  =  |x|/x,  i.e.,  sign(x)  =  l,  x>0,  sign(jc)  =  -l,  x<0. 


(30) 
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The  Fourier  transform  of  the  surface  illumination  is: 


E  ,=  E 

wk  * 


In S(k)  -  —sin{wk) 


(31) 


The  ^-dependent  radiance  distribution  just  below  the  sea  surface  is: 

8{k)  -  sin{wk)/ikn) 


ux,o,n,<p)= 

+H) 


(l  +  litT) 


p(cosYs)  + 


(On 


2(\  +  ikr) 


¥M<P) 


e'^dk  .(32) 


By  taking  appropriate  integrals  we  have  the  following  solution  for  the  upwelling  radiance 
distribution  just  under  the  sea  surface: 

L{x,Q,^i,(p)  =  -  ^(H', a:.t)]  picosy,)  +  ^[l  -  F,(w,x,T)]y/,i^,(p)^,  (33) 


F^{w,x,r)  =  - 


^  |w+jc|> 

sign(w  +  x) 

l-e'l“ 

+  signin’  -  x) 

l-e'  ^  ^ 

1  J 

K  J 

f  iH’-jrl 

signT 

!  T  1 

1 

J 

=  F;(h',X,T)  +  ;;^ 
2t 

|w- 

1 

2|t| 

(w- 

_  ^  1  -  COS<jO 


I  T  i 


\w-\-x\e 


w+x 

r 


(34) 


(35) 


Far  from  the  shadow  in  each  direction  (at  jw  -  jc|,  \w  +  x\ »  |t|)  F^—>0  and  >  0  and  Eq.  (33) 
coincides  with  the  equation  for  upwelling  radiance  (18)  for  homogeneous  illumination. 

Equations  (33)-(40)  include  functions  and  x^,  given  by  Eqs.  (20)-(21),  that  involve 
integrations  over  angular  variables  p  and  <p.  In  order  to  simpUfy  these  expressions  let  us 
substitute  phase  functions  inside  integrals  in  Eqs  (20)-(21)  by  their  transport  equivalents: 

picosy)  28  +  2(1  -  2B)3(p  -  p')S(<p  -  (p'),  (36) 

p(cosy;) -+2B+  2(\-2B)S(p'-Ps)S((p'-<Ps),  (Ps  =  0,  (37) 

B  =  0.5f^p(p)dp.  (38) 


In  this  case  after  appropriate  integrations  we  have: 
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=  2Rp(^,jU')  +  2(1  -  2B)p{cosYs), 


1 

pip,in')  =  —  J  p(cosY)d(p', 


H  f^s+m 


RADIANCE  ABOVE  THE  SEA  SURFACE 

In  order  to  calculate  radiance  distribution  above  the  sea  surface  as  a  function  of  viewing  angles,  we 
have  to  do  the  following:  a)  we  should  take  into  account  the  transmission  by  sea-air  interface  from 
below  by  multiplying  result  by  the  transmission  coefficient: 

Tjip)^l-Rl(py,  (42) 

b)  we  should  express  cosine  \p\  through  the  zenith  viewing  angle  6 : 

\p\  =  ^l-sm^e/nl-  (43) 

It  means  that  we  should  make  the  following  two  substitutions  in  Eq.  (33): 

0  0  i  t  [i-/?J(z,)1[i-r;(0)1 

El  ^  El  T,\z,)tJ ip)  ^  E,  i - iL - J ,  (44) 

cosZ^ 

and 

p^^l-sin^  d/nl-,  (45) 

here  and  R^  are  Fresnel  reflection  coefficients  of  sea  surface  for  the  light  falling  from  above 
and  below: 


KW- 


1  ( 

2  U  +  [kI^  +  P„j 


2  {n„cos0  +  pj  [cosG  +  n^pJ 
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APPROXIMATION  FOR  THE  PHASE  FUNCTION  OF  SCATTERING 


The  backscattering  coefficient  of  Ught  in  seawater  is  connected  with  the  average  cosine  of  Ught 
distribution  by  the  following  equations  (Haltrin,  1985, 1998a;  Haltrin  and  Kattawar,  1993): 


ax 

1^' 


1  +  4//' -Ji* 


(48) 


This  regression  is  derived  by  the  author  from  the  experimental  results  by  Timofeyeva  (1971): 
Ji  =  y(2.6178398  +  y(-4.6024180  +  y(9.00406(X)  +  y(-14.59994  + 


+  y(14.83909  +  y(-8.1 17954  + 1.8593222 y))))),  y= 


(49) 


An  empirical  representation  of  the  Petzold’s  experimental  angular  scattering  coefficient  p(d)  and  of 
the  scattering  phase  function  p(0) ,  where  6  is  the  scattering  angle  in  degrees,  are  represented  by 
the  following  equations  (Haltrin,  1997): 


Pi0)  =  exp 


9 


q  =  2.598+  17.748  Vh  - 16.722 f>  +  5.932 , 

=1,1 88  “  0.688  Q)q  , 

1^2  =0.1(3.07 -1.90  tOo), 

ikj  =0.01(4.58 -3.02 tUo),  ► 
=0.001(3.24 -2.25  (Uo). 

k,  =  0.000 1(0.84  -  0.6 ItUo), 

Ji+X(-i)'’^„0O 

V  n  =  l  J 


p(0)  =  — exp 
b 


(50) 

(51) 


(52) 


(53) 


The  strong  regressions  given  by  Eqs.  (50)-(53)  can  be  used  as  a  basis  for  the  empirical  model  of 
the  phase  functions  with  the  coefficients  dependent  on  the  absorption  and  scattering  coefficients. 
The  single-scattering  albedo  used  here  varies  from  0.09  to  0.96. 
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EXAMPLES  OF  COMPUTATIONS 

To  illustrate  presented  theory  and  algorithm,  several  examples  of  the  sea  shadowed  by  a  long  pier 
are  presented  here.  All  calculations  are  made  with  the  FORTRAN  program  RADPIER  listed  in 
Appendix  A.  Below  we  present  one  line  graph  (Fig.  2)  that  illustrates  the  range  of  changes  in  the 
reflected  radiance  inflicted  by  the  pier  and  six  density  graphs  (Figs. 3-8)  that  illustrate  the  two- 
dimensional  pictures  of  reflected  ra^ances. 

The  one-dimensional  graph  presented  by  Fig.  2.  shows  radiance  distribution  of  reflected  light  over 
Ox  axis  in  the  direction  perpendicular  to  the  pier.  The  values  of  inherent  optical  properties  are  taken 

as  follows;  c  =  0.5  m~\  b  =  0.2  m“';  the  solar  zenith  angle  is  taken  to  be  45°,  and  the  solar  azimuth 
angle  is  equal  to  0°  (the  solar  plane  is  orthogonal  to  the  pier  direction);  the  viewing  zenith  angle  is 
set  to  30°.  Here  we  considered  two  viewing  cases:  (1)  viewing  azimuth  angle  equal  to  15°,  and  (2) 
viewing  azimuth  angle  is  equal  to  105°.  In  Fig.  2  we  see  a  very  sharp  shadow  with  a  peculiar 
spatial  behavior  that  is  a  result  of  backscattering  from  the  sea.  Zero  values  of  radiance  in  this 
picture  are  results  of  neglecting  the  part  of  illumination  that  is  due  to  the  diffuse  hght  of  the  sky. 

Next  two  figures  (Fig.  3  and  Fig.  4)  show  two-dimensional  density  plots  of  normahzed  sea 
radiance  as  a  function  of  coordinate  x  and  azimuth  angle  of  viewing.  Zenith  angle  of  viewing  is 
30°,  sol^  zenith  angle  is  45°.  Fig.  3  displays  the  radiance  with  progressive  gray-scale  palette  and 
Fig.  4  displays  the  same  values  as  Fig.  3  only  plotted  with  the  banded  gray-scale  palette  to  show  a 
fine  structure  of  the  hght  field. 

The  next  two  figures  (Fig.  5  and  Fig.  6)  display  two-dimensional  density  plot  of  normahzed  sea 
radiance  as  a  function  of  coordinate  x  and  azimuth  angle  of  viewing.  Zenidi  angle  of  viewing  is 
35°,  solar  zenith  angle  is  45°.  Fig.  5  displays  the  radiance  with  progressive  gray-scale  palette  and 
Fig.  6  displays  the  same  values  as  Fig.  5  only  plotted  with  the  banded  gray-scale  palette  to  show 
the  fine  structure  of  the  hght  field. 


The  last  two  figures  (Fig.  7  and  Fig.  8)  are  similar  to  the  previous  figures  only  for  different  angles: 
zenith  angle  of  viewing  is  75°,  and  solar  zenith  angle  equal  to  45°.  The  difference  in  figures  7  and  8 
are  also  in  palettes:  Fig.  7  is  drawn  with  progressive  gray-scale  palette,  and  Fig.  8  is  drawn  with 
the  banded  gray-scale  palette  to  show  a  fine  structure  of  the  hght  field. 


CONCLUSIONS 

The  equations  (^)-(47)  expresses  the  upwelhng  sea  radiance  distribution  as  a  function  of  solar 
zenith  angle,  azimuth  and  zenith  angles  of  viewing,  the  width  of  the  pier  and  the  inherent  optical 
properties  of  the  water  (absorption  and  attenuation  coefficients  and  the  phase  function  of 
scattering).  Results  for  the  actual  radiance  distribution  of  hght  are  obtained  numerically  using  the  in 
situ  measured  optical  properties  of  seawater  and  are  shown  in  Figs.  3-8.  The  scattering  phase 
function  used  in  these  calculation  is  estimated  through  the  regression  equations  (50)-(53)  that  are 
derived  from  an  analysis  of  ah  fifteen  Petzold  scattering  phase  functions.  The  derived  scattering 
phase  function  depends  the  scattering  coefficient,  the  single-scattering  albedo,  and  the  scattering 
angle. 

The  main  final  result  of  this  report  consists  of  the  Eqs.  (33)-(35)  for  the  upwelhng  sea  angular 
radiance  distribution  as  a  function  of:  1)  the  solar  zenith  angle,  2)  azimuth  and  zenith  angles  of 
viewing,  3)  the  width  of  the  opaque  body,  3)  absorption  and  attenuation  coefficients,  and  4)  the 
phase  function  of  scattering. 


Ratio  of  radiance  with  the  shadow  to  the  backgound  radiance 
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Coordinate  x  orthogonal  to  the  pier,  in  m 


Coordinate  x  orthogonal  to  the  pier, 
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Fig.  3.  Two-dimensional  density  plot  of  normalized  sea  radiance  as  a  function  of  coordinate  x  and 
azimuth  angle  of  viewing.  Zenith  angle  of  viewing  is  30°,  solar  zenith  angle  is  45°. 


Coordmiftte  x 
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Fig.  4.  Two-dimensional  density  plot  of  normalized  sea  radiance  as  a  function  of  coordinate  x  and 
azimuth  angle  of  viewing.  Zenith  angle  of  viewing  is  30°,  solar  zenith  angle  is  45°.  The 
same  as  Fig.  3  only  plotted  with  the  banded  gray-scale  palette  to  show  the  fine  structure  of 
the  light  field. 


Coordinate  x  orthogorial  to  the  pier 


Distribution  of  Light  Radiance  Reflected  from  the  Sea  Illuminated  by  the  Sun  with  Shadow  1 7 


Normalized  sea  radiance 


0.00  0.25  0.50  0.75  1.00 

Normalized  sea  radiance 


Fig.  5.  Two-dimensional  density  plot  of  normaUzed  sea  radiance  as  a  function  of  coordinate  x  and 
the  azimuth  angle  of  viewing.  Zenith  angle  of  viewing  is  35°,  solar  zenith  angle  is  45°. 


Coordinate  x  orthogonal  to  the  pier,  m 
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Fig.  6.  Two-dimensional  density  plot  of  normalized  sea  radiance  as  a  function  of  coordinate  x  and 
azimuth  angle  of  viewing.  Zenith  angle  of  viewing  is  35°,  solar  zenith  angle  is  45°.  The 
same  as  Fig.  5  only  plotted  with  the  banded  gray-scale  palette  to  show  the  fine  structure  of 
the  light  field. 


Coordinate  x  orthogonal  to  the  pier,  i 
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Fig.  7.  Two-dimensional  density  plot  of  normahzed  sea  radiance  as  a  function  of  coordinate  x  and 
azimuth  angle  of  viewing.  Zenith  angle  of  viewing  is  75°,  solar  zenith  angle  is  45°. 


Coordinate  x  orthogonal  to  the  pier. 
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Fig.  8.  Two-dimensional  density  plot  of  normalized  sea  radiance  as  a  function  of  coordinate  x  and 
azimuth  angle  of  viewing.  Zenith  angle  of  viewing  is  75°,  solar  zenith  angle  is  45°.  The 
same  as  Fig.  7  only  plotted  with  the  banded  gray-scale  palette  to  show  the  fine  structure  of 
the  light  field. 
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Appendix  A 

FORTRAN  PROGRAM  RADPIER 


C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 


program  RadPier 

Created  May  5,  1997.  Last  modification;  May  9,  1997 
Written  by  Vladimir  I.  Haltrin,  NRL  Code  7331 
Ph. :  (228)  688-4528,  e-mail;  <haltrin@nrlssc . navy .mil> 

This  code  can  be  downloaded 

from  an  anonymous  ftp  site  at  <serv7330 .nrlssc .navy .mil> 


This  program  calculates  an  x-profile  of  radiance  reflected  by  the 
homogeneous  sea  illuminated  by  the  direct  sun  light  and  shadowed 
by  the  pier  of  2w  width  (in  Ox  direction)  and  infinitely  stretched 
in  the  Oy  direction.  The  depth  is  measured  along  the  Oz  axis. 

The  radiance  L  is  a  function  of  two  viewing  angles  thet  and  fi, 
and  horizontal  coordinate  x  (which  is  perpendicular  to  the  pier) . 

L  is  also  parametrically  dependent  on  the  inherent  optical 
parameters  of  water  (b,  c) ,  the  sun  zenith  angle  zang,  and  the 
halfwidth  of  the  pier  w. 


This  program  is  based  on  the  two-dimentional  in  space  (x,  z) 
theoretical  algorithm  derived  in  April  1997  by  V.I.  Haltrin 
specifically  for  the  cooperative  task  No.  5979-A7  with  the  NRL 
Code  7240  (now  Code  7340)  (P.I.  Dr.  Peter  M.  Smith). 

The  algorithm  is  based  on  the  scalar  radiative  transfer  theory. 
It  uses  Green's  function  approach  to  solve  the  problem  in 
z-direction  and  the  Fourier-transform  approach  to  solve  the 
problem  in  x-direction.  In  some  integrals  in  order  to  receive 
simple  analytical  solutions  the  phase-f unctions  are  approximated 
by  the  transport  phase  function. 

It  is  possible  to  rewrite  this  code  without  such  approximation; 
but  in  that  case  it  will  incorporate  two-dimentional  numerical 
integrations . 


implicit 

none 

integer 

i,  j  ,M,N,Nx,Na 

parameter 

(N=501,M=181) 

real*8 

x(N) ,L{N,M) , fia(M) 

real*  8 

a, c , w, zang, thet , f i , xO , dx 

real*8 

b, omeg, bp, bb,mua,mus, Rb, RO , theta, mu, xx 

real*8 

psip, f sip, angd, fourpi , phf ,mux, f is 

real*8 

FI , F2 , tau, op, LO , rLll , xi , df i 

real*8 

muav,muwa, xgor , fpsip, cosgam, pphf , ftau 

character 

tb 

data 

fourpi  712.5663706/ 

functions 
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open { 11 ,  f ile= ' radpier .in',  status^ ' old ' ) 

read{ll,*)  b  !  water  scattering  coefficient  in  1/m 

read{ll,*)  c  !  water  attenuation  coefficient  in  1/m 

read(ll,*)  w  !  pier  halfwidth  in  m 

read (11,*)  zang  !  sun  zenith  angle  in  degrees 

read (11,*)  thet  !  viewing  zenith  angle  in  degrees 

read (11,*)  Na  !  Number  of  viewing  azimuth  angles  in  degrees 

readdl,*)  Nx  !  size  of  x-array  (Nx  <  N) 

readdl,*)  xO  !  initial  x  coordinate  in  m 

read(ll,*)  dx  !  spatial  step  in  positive  Ox  direction  in  m 
close ( 11 ) 


a  =  c-b  !  absorption  coefficient 

omeg  =  b/c  !  single-scattering  albedo 

mua  -  muav(omeg)  !  average  cosine 

mus  =  muwa(zang)  !  cosine  of  the  sun  zenith  angle  in  water 
mu  =  muwa(thet)  !  cosine  of  zenith  viewing  angle  in  water 
XX  =  xgor(mua)  !  Gordon's  coefficient  x=bb/ (a+bb) 

bb  =  a*xx/(l.-xx)  I  backscattering  coefficient 

bp  =  bb/b  I  backscattering  probability 

psip  =  fpsip (bp,mu, mus)  !  auxilary  function  of  the  theory 
LO  =  mus*omeg/ f ourpi 
LO  =  L0/(mus+mu) 
op  =  0 . 5*omeg*psip 

fis  =0.  !  the  sun  plane  is  orthogonal  to  the  pier 

dfi  =  180./(Na-l)  !  increment  in  viewing  azimuth  angle 

do  j  =  1 ,  Na 

fi  =  dfi* ( j-1) 

fia(j)  =  fi  I  viewing  azimuth  angle 

mux  =  cosgam(mu,mus, f i , f is)  !  cos  of  direct  sun  rays  in  water 

angd  =  ACOSD(mux) 

phf  =  pphf (b, omeg , angd)  !  scattering  phase  function 

tau  =  ftau (c,mu, fi) 
do  i  =  1 ,  Nx 

xi  =  xO+dx* ( i-1 ) 

X ( i )  =  xi 

call  f f 12 (w, xi , tau ,  F1,F2) 

L(i,j)  =  LO* ( ( 1 . -FI ) *phf + ( 1 . -F2 ) *op)  !  Total  radiance 


end  do 


end  do 


rLll  =  l./L(l,l) 
tb  -  CHAR (9) 

open(21,  f i le= ' radpier . inf ' ,  status^ ' new ' ) 

write{21,55)  'Distribution  of  radiance  reflected  by  the  sea' 
write(21,55)  'shadowed  by  long  pier  oriented  along  axis  Oy . ' 

write  (21,55)  ’ - ' 

write(21,66)  'water  scattering  coefficient  in  1/m,  b= ' , b 
write(21,66)  'water  attenuation  coefficient  in  1/m,  c=’,c 

write(21,66)  'pier  halfwidth  in  meters,  w= ' ,  w 

write (21, 66)  'sun  zenith  angle  in  degrees,  zang= ' , zang 

write(21,66)  'viewing  zenith  angle  in  degrees,  thct=',thet 

write(21,66)  'Function  phf  =  pphf (b, omeg, angd)  =  ',phf 

write (21, 66)  'Function  op  =  0 . 5*omeg*psip  =  ' , op 

write(21,66)  'initial  x  coordinate  in  meters,  x0:=',x0 

write(21,66)  'spatial  step  in  positive  Ox  direction,  dx=',dx 
close (21 ) 


# 
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open (22 ,  f ile= ' radpier .out ' ,  status= ' new' , recl=1800 ) 
write (22,73)  Nx, tb,Na 
write(22,73)  0.,tb,0. 

write(22,75)  x(l),(tb,  x(i),  i=2,Nx) 

write (22, 77)  fia(l) , (tb, fia(j) , j=2,Na) 
do  i  =  1 ,  Nx 

write  (22, 88)  rLll*L(i,l) ,  ( tb, rLll*L ( i , j ) , j=2,Na) 
end  do 
close (22 ) 

55  format (a4 5) 

66  format (a42 , f 8 . 3 ) 

73  format (i4,al, i4) 

75  format ( f 8 . 2 , 200 (al , f 8 . 2 ) ) 

77  format(f8.2,180(al,f8.2) ) 

88  format (f8.6,180(al,f8.6) ) 

end 


real*8  function  muav (omega) 

Q  ★★***■*■★***■*•*******★★***★********★***★★★■*■****★•*■**★**★★*★*★•*■■*•*★****★★*** 

c  This  regression  is  derived  by  the  author  (V. I . Haltrin)  from  the 

c  experimental  results  by  V.  A.  Timofeyeva,  "Optical  characteristics 

c  of  turbid  media  of  the  seawater  type,"  Izv. ,  Atm.  Ocean.  Physics, 
c  Vol.7,  No. 12,  pp. 1326-1329  (English  translation:  pp. 863-865),  1971. 

implicit  none 
real *8  omega,  y 

y  =  SQRT ( 1 . - omega ) 

muav  =  y* (2. 6178398+y* (-4. 6024180+y* (9.0040600+ 

&  y* (-14. 59994+y* (14. 83909+y* (-8. 117954+1. 8593222 *y) ) ) ) ) ) 

return 

end 


Q  ★•*•**★*★★★**★★★**★★*★★*★*★*★****★★*★★****★★*★*************★***■****★★★** 

real *8  function  muwa(thet) 

Q  *****★**★**★***★**★******★******★****★★*★**★**********  +  **★*★★*★**★**** 

c  Here  a  text-book  formula  is  coded. 

implicit  none 
real *8  thet,z,nw 

data  nw  /I. 34/ 

z  =  SIND(thet) /nw 
muwa  =  SQRT ( 1 . - z  *  z ) 

return 

end 
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realms  function  xgor(inua) 

c  Formula  used  here  is  taken  from  the  book  chapter  by  V.I.  Haltrin: 

c  "Propagation  of  Light  in  a  Sea  Depth,"  Ch.  2  in:  Remote  Sensing  of 

c  the  Sea  and  the  Influence  of  the  Atmosphere  (in  Russian) ,  Moscow- 

c  Berlin-Sevastopol ,  Publ .  by  the  GDR  Academy  of  Sciences  Institute 

c  for  Space  Research,  pp.  20-62,  1985. 

Q  ★*★*★****★★***★*★★★★★★★*★★★*★★★★★★★***★★*★★★*★*★*★*★*★★*★******★*★★★★* 

implicit  none 
real *8  mua,mu2,d 

mu2  =  mua*mua 
d  =  1 . -mu2 
d  =  d*d 

d  =  d/ (1 .  *fmu2*  (4  . -mu2 )  ) 
xgor  =  d 

return 

end 


# 


real*8  function  fpsip (bp, mu, mus) 

c  Formula  used  here  is  derived  by  V.I.  Haltrin  specifically  for  the 

c  task  No.  5979-A7  for  the  NRL  code  7240  (P.I.  Peter  Smith) 

Q  ie*ic-kit**'kifk***itiiir*ir*ir*-kifkititit'k**’kiriritit**-k*ifk-k-kir*ir*ir1rir1ririririritirit'kitirir’kiritit*itirit 

implicit  none 

real *8  bp, mu, mus, Im, fp,ps,dm 


Im  =  LOG(l.+  l./mu) 

Im  =  mu*lm 
fp  =  1 . -bp 
dm  =  fp-bp 
ps  =  bp*fp*lm 

ps  =  ps+ (bp*fp*mus+dm*dm*mu) / (mus+mu) 
fpsip  =  4 . *ps 


return 

end 


real* 8  function  pphf (b, omeg, angd) 

Q  ***★★★★★*★★★★★★★★*★★★★★★★★★★★★★★★★★***★**★*★*★★*★**★*★**★★★★★****★**★★ 

c  The  equations  for  this  function  are  taken  from:  V.  I.  Haltrin, 

c  "Theoretical  and  empirical  phase- functions  for  Monte-Carlo  calcu- 

c  lations  of  light  scattering  in  sea  water."  in:  Proceedings  of  the 

c  Fourth  International  Conference  Remote  Sensing  for  Marine  and 

c  Coastal  Environments:  Technology  and  Applications,"  Vol .  I, 

c  Publication  by  Environmental  Research  Institute  of  Michigan, 

c  Ann  Arbor,  Michigan,  USA,  pp.  509-518,  1997. 

c  These  equations  approximate  all  15  Petzold  phase  functions 

c  and  relate  them  with  the  scattering  coefficient  b  and  single- 
c  scattering  albedo  omeg;  angd  is  the  scattering  angle  in  degrees. 

Q  *★*★*★★★*★★*★*★*★**★★★**★**★★★♦*★*★★★★★♦★♦★**★★★★★★*★*★****★*★★*★★*★** 
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implicit  none 

real*8  b, omeg, angd, sa, sb, q 

real*8  kl , k2 , k3 , k4 , k5 , sm, bt , pi 

data  pi  /3. 14159265/ 

sa  =  SQRT(angd) 
sb  =  SQRT(b) 

q  =  2.598+  sb*  (17.748+sb* {-16.722  +  5.932*sb) ) 

kl  =  1.188-0. 688*omeg 

k2  ^  0.1* (3.07-1.9*omeg) 

k3  =  0.01* (4.58-3 .02*omeg) 

k4  =  0.001* (3. 24-2. 25*omeg) 

k5  =  0.0001* (0.84-0. 61*omeg) 

sm  =  sa* (-kl+sa* (k2+sa* (-k3+sa* (k4-k5*sa) ) ) ) 

bt  =  EXP(q* (l.+sm) ) 

pphf  =  (4.*pi*bt)/b 

return 

end 


real*8  function  cosgam(mu,mup, f i , f ip) 
c  Here  a  text-book  formula  is  coded. 

Q  *+**★*★★★**★*★★★★+***★*★***★****★**★*****★*★****★*★**★★*★★★*★★+★★★★*** 

implicit  none 

real *8  mu, mup , f i , f ip, eg 

eg  =  COSD(fi-fip) 
eg  =  cg*SQRT (1 . -mup*mup) 
eg  =  cg*SQRT (1 . -mu*mu) 
eg  =  cg+mu*mup 
cosgam  =  eg 

return 

end 


for  the  subroutine  ffl2 


tau  =  SQRT ( 1 . -mu*mu ) 
tau  =  tau*COSD(fi) 
tau  =  tau/c 

if  (tau.eq.O)  tau  =  l.e-99 
ftau  =  tau 

return 

end 


Q  ********************************* 

real*8  function  f tau (c , mu, f i ) 

c  This  is  auxilary  function 

implicit  none 
real*8  c,mu,fi,tau 
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subroutine  f f 12 (w, x, tau,  fl,f2) 

c  Formulae  used  here  is  derived  by  V.I.  Haltrin  specifically  for 

c  the  task  No.  5979'-A7  for  the  NRL  Code  7240  (P.I.  Pit  Smith) 

implicit  none 

real *8  w, x, tau, f 1 , f 2 ,  wl , w2 , si , s2 , st , el , e2 , al , a2 , at , dl , d2 


wl 

=  w-x 

w2 

=  w+x 

si 

=  1. 

if 

(wl . le. 0 . ) 

si  = 

-1 

s2 

=  1. 

if 

(w2.1e.O. ) 

s2  = 

-1 

st 

=  1. 

if 

(tau. le. 0 . ) 

St  = 

-1 

el  =  ABS(wl/tau) 
el  =  EXP (-el) 
e2  =  ABS(w2/tau) 
e2  =  EXP(-e2) 

dl  =  sl*(l.-el)+s2*(l-e2)+st*(e2-el) 

fl  =  0.5*dl 

al  =  ABS(wl) 

a2  =  ABS(w2) 

dl  =  (al*el-a2*e2 ) /tau 

at  =  ABS(tau) 

d2  =  (wl*el+w2 *e2 ) /at 

f2  =  fl  +  0.5*(dl-d2) 

return 

end 


A  sample  of  an  input  file  'radpier.in'  (starts  just  after  the  dotted  line): 


0.2b  =  water  scattering  coefficient  in  1/m 

0.5c  =  water  attenuation  coefficient  in  1/m 

1.5  w  =  pier  halfwidth  in  meters 

45.  zang =  sun  zenith  angle  in  degrees 
25.  thet =  viewing  zenith  angle  in  degrees 


181  Na 
201  NX 
-10. 
0.1  dx 


xO 


=  number  of  viewing  azimuth  angles  in  degrees  (0-7C) 
=  size  of  x-array  (Nx  <  N) 

=  initial  x  coordinate  in  meters 
=  spatial  step  in  the  positive  Ox  direction  in  m 
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